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Abstract. The thin-layer Navier-Stokes equations are coupled with a 
zonal scheme (or domain-decomposition method) to develop the Transonic Navier- 
Stokes (TNS) wing-alone code. TNS has a total of 4 zones and is extended to a 
total of 16 zones for the wing-fuselage version of the code. Results are 
computed on the Cray X-MP-48 and compared with experimental data. 

Key words, zonal methods, applied aerodynamics, computational fluid dynamics, 
transonic, viscous. 

Running head: Zonal-Method Simulation of Transonic Flows 

1. Introduction. With the improvement of numerical algorithms for the 
solution of the three-dimensional Euler/Navier-Stokes equations and the recent 
advancements in computer capabilities, previously unchallenged problems in 
computational fluid dynamics (CFD) are now being attempted. Some of the more 
recent transonic applications involving more sophistication in the geometries 
have appeared; they include Deiwert and Rothmund [9], Deiwert et al. [10], 
Fujii and Kutler [15], and Mansour [24]. However, even these solutions were 
conducted on coarse grids and required large amounts of cpu time, thereby 
precluding their use for more complicated geometries. 
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The next step in the advancement of CFD is in the simulation of transonic 
flow over complex geometries. However, generating a single grid in three 
dimensions for this simulation is one of the biggest pacing items in CFD. 
Currently, numerically generated grids in curvilinear coordinate systems are 
commonly used around arbitrary geometries. The implementation of the boundary 
conditions is easier on body-fitted systems. As the geometry becomes more 
complex, such as wing-body-strake or wing-canard combinations, the generation 
of a single grid becomes a very difficult task. Add to this the requirement 
of appropriate clustering of the grid points at all no-slip surfaces, and in 
regions of high gradient flow, and this task becomes nearly impossible. 

To help alleviate this problem, zonal approaches have become increasingly 
popular. In the zonal approach, the flow field is partitioned into distinct 
"zones," each of which is solved independently. Using this approach, a coarse 
base grid can be generated about the geometry in question. Then the grid for 
any subdomain can be generated easily. Refined zones for regions of high 
gradients such as shear layers, shock waves, jets, wakes, and vortices can 
also be easily generated. As more components are added to the geometry, more 
zones can be added to capture the pertinent flow about them. 

This zonal approach is not novel, having been used implicitly in the 
introduction of boundary-layer theory by Prandtl [25]. Here the flow field is 
divided into an inner viscous region where the boundary -layer equations are 
solved and matched with an outer inviscid solution. Thus, using this zonal 
approach, the boundary-layer equations could be solved via a marching tech- 
nique while different relaxation methods are applied to the outer inviscid 
flow. Most noncomplex general viscous flows can be simulated without using 
the more consuming Navier-Stokes equations. 
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The zonal approach has a number of advantages: (i) the difficulties in 

generating three-dimensional grids for different types of complex configura- 
tions can be reduced with the use of zonal methods; (ii) different types of 
grid topologies can be implemented where appropriate in order for the grids to 
be mesh-efficient (that is, more points on the body surface where accuracy is 
desired, and fewer points in the flow field); (iii) different equations sets, 
with different algorithms, can be implemented in the various zones for compu- 
tational efficiency; and (iv) information for only one zone need reside in the 
computer core at any given time, thereby relaxing memory limitations. Even in 
turbulence modeling, the zonal approach can be applied as described by Kline 
et al. [21]. In short, there is a wide variety of situations in which zonal 
methods can be applied with substantial gain. 

Figure 1 illustrates the two major types of zonal interfaces, patched and 
overlapped. Earlier work in the development of a zonal -boundary scheme, for a 
system of hyperbolic equations includes that of Cambier et al. [6]. Further 
work on patched grids using the Euler equations was done by Rai [30], 

Hessenius and Pulliam [16], and Hessenius and Rai [17]. The patched grid 
concept was also used for the potential and full -potential equations by Lee 
et al. [22] and Yu [35]. The approach of overlapped grids was used with the 
Euler equations by Atta [1] and Atta and Vadyak [2]. Steger et al. [33] give 
results obtained on overlaid grids in conjunction with the stream- function 
approach. The current work uses the patched grid zoning method in a grid- 
refinement mode. Other work using this grid refinement mode has been done by 
Eriksson [11] and Baker et al. [3]. It is also possible with the zonal 
approach to solve different types of equation sets in the different zones [8]. 

This paper discusses the computational approach in which the fast- 
convergent, Pulliam-Chaussee [28] diagonal algorithm is coupled with a zonal 
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approach. The new approach permits relatively inexpensive fine-grid solutions 
to be made of the Euler/Navier-Stokes equations, which is especially important 
for flows with shock/boundary-layer interaction. Validation of the code by 
comparing numerical solutions with existing experimental data will also be 
presented. 


2. Governing equations. The equations solved in this study are the 
Reynolds-averaged Navier-Stokes equations written in strong conservation-law 
form. These equations are simplified by using the standard thin-layer approx- 
imation for the viscous terms. The thin-layer Navier-Stokes equations written 
in generalized curvilinear coordinates are 
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with 
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Pressure is related to the conservative flow variables Q by the equa- 
tion of state 

p = (y - 1 ) Je - ^ p(u 2 + v 2 + w 2 )l . 


The Beam-Warming algorithm [5] is used to solve the governing equations, 
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where A, B, C, and M are the Jacobian matrices 3E/3Q, 3F/3Q, 3G/3Q, and 
3S/3Q, respectively. Note that M, which is derived from S, contains 
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derivatives in n. For h = (1/2) or 1, the integration is trapezoidal 
(second-order) or Euler implicit (first-order), respectively, in time. These 
equations are central-space-differenced and implicitly advanced in time. To 
maintain stability of the algorithm (because of the central-difference scheme 
used), an explicit fourth-order artificial-dissipation term is added to the 
flux calculations and an implicit second-order dissipation term is added to 
each of the block tridiagonals. The first working code used the above 
algorithm. 

For steady-state computations or first-order time- integrations, a diag- 
onal form of (2) can be used. In this case, the left and right eigenvector 

A /V A 

matrices of A, B, and C are used to diagonalize the one-dimensional opera- 
tors. The diagonal algorithm in three dimensions has the form 

(3) T (I + hi A )N( I + hi A )P(I + hi A )T -1 AQ n = R n . 

e, ( ( nn ccc 

For a complete derivation of the diagonal algorithm (as well as definitions 
for N, P, etc.) see [28]. 

The main advantage of this form is the simplification of the matrix 
inversions from block-tridiagonal inversions to scalar-tridiagonal inver- 
sions. This simplification reduces the computational work by about 30^ 

[28]. Also, the new scalar form for the inversion process allows the use of 
scalar-pentadiagonal solvers so that the added fourth-order explicit artifi- 
cial dissipation can be properly linearized and be made fully implicit. This 
form enhances stability and convergence rates [26,29]. To further enhance the 
convergence rate, a space-varying At has been used. It is given by the 

1 /p 

formula At = AW[1 + (J) ], where At Q , as used here, is simply a constant 

used to decrease or increase At. For most cases, the default value of At Q 
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was 5. In viscous calculations the diagonal algorithm uses an explicit treat- 
ment of the viscous terms. The turbulence model used for all cases is the 
Baldw in-Lomax algebraic model [4]. 

3. Zonal approach. To initiate the zonal approach, a coarse block 
(zone 1) is first generated about the configuration. The coarse block can be 
generated either iteratively, elliptic [32], or through a marching scheme, 
parabolic [12]. Both procedures have the capabilities of spacing and orthog- 
onality at the inner and outer boundaries. The topology of the grid is H-type 
in both the spanwise and chordwise directions. 

To generate the finer zones near the wing, a small zone of points about 
the wing is removed from grid 1 . The space left open by the removal of points 
from zone 1 is then occupied by the finer grid (zone 2). Zone 2 is generated 
by putting twice as many points in every spatial direction relative to 
zone 1 . This task is accomplished by cubic-spline interpolation of the 
coarse-grid points to the fine-grid points. To generate the viscous grids, a 
small zone of points is again removed about the wing from zone 2. Zones 3 
and 4 now occupy the area left vacant by the removal of points from zone 2. 
Zone 3 occupies the area above the wing and includes the upper surface of the 
wing, and zone 4 occupies the area below the wing and includes the lower 
surface of the wing. Zones 3 and 4 retain the same number of points in the 
streamwise and spanwise direction as does zone 2; however, points are further 
clustered in the normal direction to capture viscous effects. All zones 
overlap at the zonal boundaries, usually by one or two grid planes. 

Figure 2 shows a typical grid with outer boundary positions specified so 
as to coincide with the position of the wind-tunnel walls from the NASA Ames 
High Reynolds Number Channel I [231. The grid is plotted in perspective so 
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that detail on the upper and lower wind-tunnel-wall surfaces, the inflow and 
outflow planes, and the wing-symmetry plane are all visible. This grid, which 
is generated directly by the parabolic grid-generation approach, becomes the 
outer, coarse-grid zone (grid 1). The grid detail near the wing/symmetry- 
plane juncture has been removed. The blowup of the grid in this region 
(Fig. 3) shows the detail of the grid zones 2, 3, and 4. The wing geometry 
used in this case is composed of NACA 0012 cross sections, has a taper ratio 
of 1.0, a leading edge sweep of 20°, an aspect ratio of 3.0, and is rigged in 
the wind-tunnel-wall grid at an angle of attack of 2°. This wing does not 
have any twist or dihedral. Note that the grid immediately adjacent to the 
wing surface (grid zones 3 and 4) is highly clustered in the normal direction 
and is, therefore, appropriate for a Navier-Stokes flow solver. Also note 
that the Navier-Stokes grid expands in thickness from the leading edge toward 
the trailing edge so as to better capture the growing boundary layer. 

Figure 4 illustrates the wing surface grid and a chordwise slice of 
zones 1, 2, and 3 at the symmetry plane (y = 0). Note the doubling of points 
in zone 2 (in chordwise and normal directions) relative to zone 1. The dou- 
bling is also done in the spanwise directions, although this is not shown 
here. There can also be seen a one-to-one correspondence of points in the 
chordwise direction (also in the spanwise direction) between zones 2 and 3, 
and the clustering of grid cells in zone 3. Also highlighted are the wing 
surface, trailing-edge mesh, and wing-tip region. Zone 4 is not shown, but it 
is essentially of the same structure as zone 3. 

For clarification, Fig. 5 shows a generic form of the overlapping proce- 
dure between two inviscid zones (this generic overlap procedure is the same 
between all zones, inviscid or viscous). Boundary conditions are applied 
explicitly when using this overlapping procedure. That is, boundary 
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conditions for the zone 2, the £ = A plane, are obtained by interpolating 
data from the interior points of zone 1. And, conversely, boundary conditions 
for the zone 1, the 5 = B plane, are obtained from the interior data from 
zone 2. The boundary conditions for a single surface are obtained with just a 
series of one-dimensional interpolations. The interpolation process is auto- 
mated to the extent that only the two planes involved in the interpolation 
need to be defined, the base and target planes. Two interpolation schemes are 
coded: one is a cubic spline and the other is linear. The cubic spline is 
best in smooth regions of the flow, and the linear is best in the nonsmooth 
flow regions. The linear interpolation routine was used for the results pre- 
sented herein. Figure 6 illustrates the possible arrangement when more zones 
than one are used to update the boundary conditions of zone 1. Finally in 
Fig. 4, zones 1 and 2 are solved using the Euler equations (inviscid zones), 
and the Navier-Stokes equations are used in the viscous zones (zones 3 and 4). 

4. Discussion of results. 

4.1. Transonic Navier-Stokes (TNS) wing-alone code. The first case 
tested consisted of a NACA 0012 wing, subject to the following flow condi- 
tions: = 0.826, a=2°, and a Reynolds number, based on the chord, of 

8 million. This case, which is moderately difficult, involves a strong shock 
that extends from zone 3 (viscous upper-surface zone) into zones 2 and 1 
(inviscid zones). With the zonal approach, the test case was run on a 
150,000-point mesh, which for three dimensions is a relatively fine grid. 

For this case the wind-tunnel-wall effects are very significant. This 
effect can be seen in Fig. 7 where the pressure coefficient distributions from 
the TNS code, with and without the walls modeled, are compared with those of 
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experiment [23], The shift in shock position caused by the tunnel walls is 
obvious. The shock position for the case with walls is in good agreement with 
the experimental shock position. 

General agreement between the wind-tunnel-wall case and experiment is 
better inboard of mid-semispan than it is outboard. In particular, the com- 
puted upper-surface shock strength at 2y/b =0.78 is larger than that of the 
experiment. This larger strength is caused by a large boundary-layer separa- 
tion in the experimental results at this semispan location, which is not 
accurately reproduced by the computed results. A good picture of this situa- 
tion is given in Figs. 8 and 9, which show a set of computed particle paths 
(Fig. 8) and an oil-flow photograph taken from the experiment (Fig. 9). The 
experimental separation is about twice as large as the computed separation. 

The spanwise extent of the experimental separation is reasonably predicted by 
the computation, but the streamwise extent is under predicted. Some reasons 
for this discrepancy are coarse-grid and turbulent-model effects, as well as 
the sensitivity of the flow pattern to changes in the free-stream Mach 
number. Despite the difference in the size of the separation zone, the over- 
all comparison is quite encouraging. 

Figure 10 shows Mach number contours plotted in a wing cross-sectional 
plane at a semispan station of 2y/b = 0.66; the zonal boundaries are high- 
lighted. Note the smoothness with which the shock wave crosses the zonal 
interface boundary. This demonstrates that this particular interface, between 
zones 2 and 3, is implemented in a conservative manner. (It should be noted, 
however, that not all interfaces are conservative.) Generally, most of the 
other contours cross the zonal interface boundaries in a smooth and continuous 
way. Downstream in the wake where the fine viscous grid interfaces with a 
relatively coarse inviscid grid (see Fig. 3), the wake abruptly stops. This 
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is due to the difference in grid refinement (in going from the wing grid to 
the fine, but not clustered, grid 2). The power of the zonal scheme allows 
the addition of more zones in those regions where large gradient flow is 
observed. In the wing-fuselage case, a special zone is implemented downstream 
of the wing to capture the important viscous wake effects. 

The convergence rate of the diagonal version versus the block ADI version 
is illustrated in Fig. 11. The time-step used in the nondiagonal version 
was At = 0.004, which was the largest time-step possible while still main- 
taining stability of the code. The diagonal version used a variable time-step 
(as described previously). Even though the comparison is between a fixed 
time-step solution and a variable time-step solution, the main speedup in the 
diagonal algorithm is not in the variable time-stepping procedure, but in the 
proper linearization of the dissipation terms [26]. The slow rate of conver- 
gence in the nondiagonal version seems to occur in the outer inviscid zones. 
The residual in the viscous zones in the first thousand iterations drops 
fairly fast, then begins to flatten out. In 5,000 iterations, all zones have 
dropped about two orders of magnitude in the L2 norm of the residual. In 
contrast, the convergence rate of the diagonal version drops rapidly in all 
the zones. A three-order-of-magnitude drop in the L2 norm occurs in about 
400-500 iterations. This convergence rate (coupled with the decrease in 
arithmetic operation count caused by the diagonal algorithm) increases the 
speed with which solutions are obtained by a factor of 40. The faster conver- 
gence rate, as stated before, is due to the proper linearization of the 
fourth-order-explicit dissipation operator. This rate was possible in the 
nondiagonal version, but would involve inverting block pentadiagonals, which 
would substantially increase the computational cost. More time was required 
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for the same case when wind-tunnel walls were modeled, since a stronger shock 
occurs. 

We can also look at the development of lift and the number of supersonic 
points (NSP) to study the convergence characteristics. Figure 12 shows the 
evolution where the x-axis is the number of iterations; the y-axis is left 
unlabeled (since the actual values were not of significance, but the overall 
trend in the approach to the steady state is important). The NSP overshoots 
the final converged solutions at about 100 iterations, but then quickly 
approaches the converged solution. At about 200 iterations, it is within 1 % 
of the converged solution, and at 300 iterations it is within The lift 

also overshoots the final converged value, but at 200 iterations it is within 
about of the final solution, and at 300 iterations it is within 1%. If, 
instead of a three-order-of-magnitude drop in residual, convergence is based 
on 9511 of the converged lift, then a solution can be generated in about 
200 iterations, or about 18 min of cpu time for this fine-grid calculation. 

The next case presented consists of a massive, shock- induced, boundary- 
layer separation. This case was computed to ascertain the degree of robust- 
ness of the present algorithm and, in particular, the ability of the present 
zonal interface scheme to cope with large flow gradients. The geometry used 
is the same as that of the last case. The free-stream Mach number and angle 
of attack have been arbitrarily chosen to be 0.9 and 5°, respectively. (Jtili 
zation of the wind-tunnel-wall boundaries produced a "choked" solution with a 
shock wave spanning the tunnel. After several hundred iterations and a moder 
ately converged calculation, the solution diverged as expected. This result 
was a consequence of the "fixed" upstream boundary conditions forcing more 
mass flow through the tunnel than the choked condition would allow. 
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The solution was repeated with free-air boundary conditions, and then 
convergence was easily achieved. Computed particle paths on the upper wing 
surface are displayed in Fig. 13. Note that approximately half of the upper 
wing surface is separated. The solution contains several interesting fea- 
tures, including- a separation saddle point, a focus, a reattachment saddle 
point, and a node. Note also that this computed solution has a stable 
sequence of critical points on the separation line; that is, a node followed 
by a saddle point and then a focus. This case required about 2.8 hr of cpu 
time . 

Two different perspective views of the three-dimensional particle paths 
are shown in Fig. 14. Figure 14a shows a view from outboard of the wingtip, 
and Fig. 14b shows a view from behind and above the wing. The height and 
three dimensionality of the separation zone are apparent in these figures. 

The dashed particle paths move along the wing surface until the separation 
line is encountered; then they are deflected up and over the separation 
bubble, with a few of the dashed paths captured by the primary swirling flow 
at the center of the wing. The solid particle paths are more intimately 
involved with the two swirling pockets of flow, and they essentially define 
these regions. 

The position of the separation region relative to the zonal interface 
boundary is best displayed by plotting particle paths constrained to lie in 
spanwise cross-sectional planes. Two such plots are displayed in Fig. 15. 
Figure 15a shows cross-sectional particle paths for a semispan station of 
2y/b = 0.66. The separation region is large and easily extends above the 
zonal boundary from the Navier-Stokes region into the Euler region. Neverthe- 
less, the solution looks qualitatively reasonable. An enlargement of the 
separated portion of the solution is shown in Fig. 15b. From this figure it 
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can be seen that the particle paths pass smoothly across the interface bound- 
ary with no function or slope discontinuities. Thus, the primary objective of 
doing this calculation was achieved. Despite the existence of a strong shear 
gradient across the explicitly updated interface boundary, the present 
approach is capable of predicting a stable solution that is reasonably free 
from interface boundary influence. 

Figure 16 shows Mach number contours for this difficult transonic case at 
a semispan station of 2y/b =0.33. A very strong shock can be seen on the 
upper surface of the wing, as well as its smooth transition to the other 
zones. The wake region is larger because of the shock induced separation 
region. 

For more details on the zonal procedure, convergence studies, transonic 
results with different wings, and data management structure, the reader is 
referred to [13,19,20]. 

The zonal grid topology described in the previous section utilizes 
Cartesian-like grids to simplify the zonal interfacing and to maintain flow 
conservation at shocks. One drawback to this approach, however, is the 
resulting Hd-mesh singularity at the wing leading edge. This occurs at the 
interface between blocks 3 and 4, as shown in Fig. 17. It is important to be 
able to properly treat severe coordinate singularities because they naturally 
arise in realistic aircraft configurations. 

Some of the early computations with the TNS code were performed at low to 
moderate angles of attack and used simple central differencing of the met- 
rics. For these relatively simple test cases, there appeared to be no diffi- 
culty at the wing leading edge. In order to control metric truncation errors 
and insure uniform flow as an exact solution of the finite-difference equa- 
tions, the free-stream residual was subtracted from the right side of (1). 
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However, for = 0.5 and a - 10°, large oscillations in the flow variables 
developed at the H-mesh singularity, as indicated by the wing root section 
Mach number contours given in Fig. 18. It was not possible to obtain a fully 
converged solution. A nominal two-order-of-magnitude drop in the norm of 

the residual was possible only after adding large amounts of numerical 
dissipation near the leading edge. This was clarly unacceptable and would 
prevent high-angle-of-attack simulations. This problem was fixed by using 
free-stream preserving metrics as described by Pulliam and Steger [27]. 

In their work, Pulliam and Steger used special numerical metrics that 
insured uniform flow was: (1) an exact solution of the finite-difference 

equation; or (2) the free-stream subtraction described above. In the cases 
treated in [27], Pulliam and Steger found no appreciable difference between 
the two solutions or in their convergence rates, and therefore preferred the 
latter approach. However, their grids were relatively smooth and did not 
possess severe coordinate singularities such as are present in the H-type. 

Once the free-stream preserving metrics were implemented in the TNS code, 
there was no difficulty in rapidly converging to a steady state, and no addi- 
tional dissipation was necessary at the leading edge. (For more details into 
the theory of the free-stream preserving metrics, see [7].) The resulting, 
improved, Mach-number contours are shown in Figs. 19 and 20. There is a 
nominal amount of distortion at the interface of blocks 3 and 4 because the 
flow variables are obtained there using simple averages instead of from the 
governing equations. The Lj norm of the residuals for all four blocks are 
shown in Fig. 21. A three-order-of-magnitude drop in the residuals of all 
four blocks was obtained in 700 iterations. This criterion for convergence is 
usually sufficient for plottable accuracy. 
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The first set of solutions was obtained at a free-stream Mach number 


= 0.5. A comparison between TNS lift and drag coefficients and those 
obtained with a method based on the full-potential equations [18] is shown in 
Figs. 22 and 23. Unfortunately, good-quality force, moment, and surface- 
pressure data, up through maximum lift, are lacking in the literature. 

Although force and moment data provide a way to assess a code's global accu- 
racy, extensive surface-pressure data are necessary for validating the details 
of the flow simulation. These flow simulations indicate that the flow is 
subcritical owing to the effects of three-dimensional relief and wing sweep. 
Sonic flow is achieved only at maximum lift (a = 13-5°). The lift coefficient 
obtained with the TNS code is in good agreement with the full-potential result 
in the low-angle-of -attack range, but differs significantly in the high-angle- 
of-attack range. The TNS code predicts maximum lift at a = 13.5°, whereas 
the full-potential lift coefficient continues in a linear fashion. The drag 

coefficient exhibits the usual quadratic variation with angle of attack 
2 2 

(C n - C. ~ a ) together with the large drag rise at stall. A converged solu- 
tion usually requires 700 iterations, or 55 min of Cray X-MP time. 

Particle trajectories for a = 15° are shown in Figs. 24 and 25. 

Figure 24 is a perspective view from above the wing and looking downstream 
toward the wing leading edge. The particles are released along the wing 
leading edge and wing tip. This massively separated steady flow exhibits 
a- induced separation, that is, separation caused solely by angle of attack, as 
apposed to shock-induced separation. The vortical structure of the separated 
region is evident. The particle trajectories emanating from the wing tip 
indicate a wing-tip vortex. Figure 25 is an end view of the wing looking 
inboard from the wing tip. The region of separation extends across zonal 
boundaries in a smooth manner. This solution seems to be on the verge of 
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going unsteady, and it required twice as many iterations to achieve conver- 
gence. The spatially varying time-step tends to inhibit unsteadiness. When a 
constant time-step was used, a cyclic variation in the residuals was noted. 

All of these solutions were obtained with little difficulty, provided 
consistent metrics were used. Without these special free-stream-preserving 
metrics, solutions at angles of attack greater than 5° would not converge. 
These high-angle-of-attack solutions also demonstrate the robustness of the 
zonal approach. 

Transonic wing solutions were also obtained for = 0.8. The TNS 
C L vs a vaiation is compared with the previous subcritical case in Fig. 26. 

In this transonic flow, maximum lift occurs at a much lower angle of attack 
(a = 6°), a result of shock-induced separation. These transonic solutions 
were relatively easy to obtain, requiring 50 min per solution. 

The wing-tip vortex for the maximum lift condition, a = 6°, is shown in 
Fig. 27. A very interesting simulated oil-flow pattern on the upper surface 
of the wing is shown in Fig. 28. These are particle trajectories constrained 
to the next coordinate surface above the wing (because of the no-slip condi- 
tion on the wing). Notice the saddle and nodal point singularities on the 
upper surface of the wing. These critical points form a stable topological 
configuration as defined in (34]. There is a major separation line extending 
over most of the wing span, and it is followed by a reattachment line (shown 
by a dashed line) a short distance downstream. A second small separation 
region is also evident near the trailing edge in the vicinity of the wing tip. 

The extent of separation in these subcritical and transonic cases tends 
to be underpredicted. This was also evident for the initial test case in 
which comparisons were made between simulated and experimental surface oil- 
flow patterns. There are several reasons for this discrepancy. First, better 
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grid resolution is probably required to improve the accuracy of the computa- 
tion. Grid-refinement studies using a million-point wing grid on the Ames 
Cray 2 are in progress. Second, an improved numerical dissipation model is 
required. Blended fourth-order and second-order smoothing is used in (1) as 
described by Pulliam [26]. This blended smoothing has very good shock- 
capturing chearacteristics on a grid suitable for the Euler equations. How- 
ever, the dissipation coefficient varies as 0(1/Az) and, on fine-spaced grids 
needed for viscous computations the numerical dissipation can be as large as 
the physical dissipation. Finally, an improved and efficient turbulence model 
is required that can adequately model three-dimensional shock-induced separa- 
tion. The first two points will be easily achievable in the near future, but 
the latter still remains uncertain. 

4.2 TNS wing-fuselage code. As previously mentioned, one of the advan- 
tages of the zonal method is the ability to create grids, about complicated 
aircraft, with sufficient clustering of points on all no-slip surfaces. 

Figure 29 shows body-conforming zonal grids of a fighter aircraft configura- 
tion in the physical and computational spaces. Figure 30 illustrates the base 
grid that is used by the zoner code to create the different zones about this 
modified F-16A aircraft. This grid was generated via an elliptic grid- 
generating technique [32]. From this figure the polar type of topology about 
the fuselage can be seen. This polar grid collapses into a singular line in 
front of the fuselage nose but, because of the averaging about this singular- 
ity, there is no problem in the flow solver. A chordwise slice through the 
wing would reveal an H-mesh-type of local topology. The modified F-16A can be 
seen in Fig. 31* illustrating the exclusion of the tail assembly. From 
Figs. 31a and 31b, it can be noted that the forebody, canopy, leading-edge 
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strake, wing, and shelf regions are unmodified and are of the exact geometry 
of the F-16A. The inlet is faired over. Figure 32 shows a planform view of 
the base grid. Generating even a base grid with no clustering normal to the 
surfaces is difficult. Because of the tapering of the wing (the horizontal 
tail is not computed here), the grid lines tend to squeeze together at the 
leading edge near the tip region, then suddenly stretch out as they continue 
to the outer boundary. 

The base grid is fed into a zonal routine, which then subsequently 
creates 16 zones (by subdividing the original base grid). Through certain 
parameters in the zonal routine, the zones created are essentially of three 
types: (i) inviscid zones, (ii) viscous zones with clustering on one face of 

the zone for wing or fuselage surfaces, and (iii) viscous zones with cluster- 
ing on two adjacent faces of the zone for the wing-fuselage juncture. 

Figure 33 shows a schematic view of the 18 zones created for the F-16A. 

Zones 1 and 4 are on the fuselage and are constructed with sufficient grid 
resolution in the normal direction. Zones 8 and 9 (zone 9 is not visible) are 
on the upper and lower wing, respectively, and also have sufficient grid 
clustering. Zones 2 and 3 (zone 3 is not visible) are special, in that they 
have clustering both normal to the wing and normal to the fuselage surfaces. 
This is in order to capture viscous effects caused by the wing-fuselage junc- 
ture. Zones 12 and 13 are not considered to be no-slip surfaces; however, 
they are clustered across the wing planform plane, in order to capture wake 
(zone 12) and tip (zone 13) effects. Zone 6 is also clustered both normal to 
the fuselage and across the wing planform plane. (Viscous effects and appro- 
priate turbulence model effects are implemented for these zones and solved 
with the thin-layer Navier-Stokes equations. All other zones are solved with 
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the Euler equations.) Results will be shown only for a 16-zone version of the 
modified F-16A. 

Figure 34 shows some of the different wing-fuselage mesh topologies. The 
topology chosen was that of Fig. 34a. It allowed for clustering in the wing- 
fuselage juinction, and made it easier to implement the turbulence model. 
Figure 35 shows the clustering in the viscous zones. The viscous zone is 
based on the first three inviscid points of the body, NJC = 1-3. Then a 
single parameter, per block, determines the amount of clustering to be imple- 
mented. A different interpolation routine from that of the wing-alone code is 
required. Instead of simple one-dimensional interpolation [19], two- 
dimensional interpolation is now required because of the interfacing of sur- 
faces. This can be seen in Fig. 36, where surface ABCD from a coarse zone, 
may be used to interpolate the surface onto a fine grid. 

Pressure contours are displayed in Fig. 37 for flow conditions of 
M^ = 0.9, a = 1.69°, and a Reynolds number based on the chord of 4.5 mil- 
lion. The coalescing of contours indicates regions of high gradient flow. 
Stagnation flow occurs at the fuselage nose and at the forward base of the 
canopy. The flow then accelerates over the top of the canopy. At the leading 
edge of the wing, the flow accelerates and smoothly decelerates as it 
approaches the trailing edge. Notice how the pressure contours continue from 
the upper wing zone into the fuselage zone. Similarly, the fuselage has two 
zones, the boundary occurring approximately at the base of the back of the 
canopy. Pressure contours continue smoothly between these two zones. 

Mach contours are plotted in Fig. 38. These contours were calculated on 
planes that were about 15 grid points up from the surfaces, so as to be out of 
the boundary-layer region. The peak Mach number is about 1.2, and the con- 
tours on the wing indicate a weak shock, if any, occurring. 
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These calculations required about 2,500 iterations on a total grid size 
of about 300,000 points. This required about 10 hr of cpu time on the Cray 
X-MP-48. For results from other transonic calculations, as well as at higher 
angles of attack, the reader is referred to [14,31]. 

5. Conclusions. A fast diagonal algorithm has been successfully imple- 
mented within the framework of a zonal approach. Results indicate that the 
modified code (in obtaining a solution for a moderately difficult case) still 
maintains its fast convergence characteristics. This improvement is demon- 
strated by producing a three-dimensional, fine-grid Euler/Navier-Stokes wing 
solution in 45 min on the Cray X-MP, which is faster by a factor of 40 than 
the original code. Computed pressures and oil-flow patterns compare favorably 
with those of experiment. The wing code has been extended to a 1 6-block 
version for the modified F-16A geometry. The zonal procedure allows appro- 
priate grid clustering normal to all no-slip surfaces, and pressure contours 
continue smoothly across the different zones. The calculations require about 
10 hr of cpu time for 300,000 grid points. 
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FIGURE CAPTIONS 


FIG. 1. Zonal interface alternatives. 

FIG. 2. Global finite-difference grid showing wing-tunnel walls. 

FIG. 3. Expanded view of embedded grid near wing. 

FIG. 4. Perspective view of embedded grid with upper symmetry plane 
(y = 0.0) and wing surface highlighted: X EE = 20°, AR = 3.0, TR = 1.0. 

FIG. 5. Grid zone interface procedure: (a) two-zone grid arrangement 

showing overlap, (b) grid point detail in the overlap region (5 = A plane). 

FIG. 6. Overlap between zones 2 and 3 removed for clarity. 

FIG. 7. Pressure coefficient comparisons: NACA 0012 airfoil section, 

AR = 3.0, A le = 20°, TR = 1.0, = 0.826, a = 2° , R e = 8 x 10 6 . 

FIG. 8. Computed particle paths on the upper wing surface: NACA 0012 

airfoil section, AR = 3.0, Agg = 20°, TR = 1.0, = 0.826, a = 2°, 

R = 8 x 10 6 . 
e 

FIG. 9. Oil-flow pattern on upper wing surface: NACA 0012 airfoil 

section, AR = 3.0, A LE = 20°, TR = 1.0, = 0.826, o = 2° , R e = 8 x 10 6 (from 

[23]). 

FIG. 10. Cross-sectional Mach-number contours: NACA 0012 airfoil sec- 
tion, AR = 3.0, A le = 20°, TR = 1.0, = 0.826, o = 2° , R g = 8 x 10 6 . 

FIG. 11. Convergence rate comparison. 

FIG. 12. Development of lift and number of supersonic points (NSP). 

FIG. 13. Computed particle paths on the upper wing surface: NACA 0012 
airfoil section, AR = 3.0, a le = 20°, TR = 1.0, = 0.9, a = 5° , 

R e = 8 x 10 6 . 
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FIG. 14. Computed three-dimensional particle paths on the upper wing 

surface: NACA 0012 airfoil section, AR = 3-0, A LE = 20°, TR = 1.0, = 0.9, 

a = 5°, R. = 8 x 10^: (a) view from behind and above wing; (b) view from 

V 

outboard of wing tip. 

FIG. 15. Copmputed cross-sectional particle paths: NACA 0012 airfoil 
section, AR = 3.0, A L£ = 20°, TR = 1.0, M m = 0.9, a = 5°, R e = 8 x 10 6 : 

(a) 2y/b = 0.66; (b) 2y/b = 0.66, expanded view of separated flow region. 

FIG. 16. Cross-sectional Mach-number contours: NACA 0012 airfoil sec- 
tion, AR = 3.0, A le = 20°, TR = 1.0, = 0.9, a = 5°, R g = 8 x 10 6 . 

FIG. 17. H-mesh singularity at wing leading edge (blocks 3 and 4): 

NACA 0012, M, = 0.5, a = 10.0°, R g = 8 x 10 6 . 

FIG. 18. Symmetry plane Mach-number contours with central-difference 

metrics: NACA 0012, = 0.5, o = 10.0°, R g = 8 x 10 6 . 

FIG. 19. Symmetry plane Mach-number contours with new metrics: 

NACA 0012, M,,, = 0.5, a = 10.0°, R e = 8 x 10 6 . 

FIG. 20. Symmetry plane Mach-number contours with new metrics: 

NACA 0012, M,,, = 0.5, a = 10.0°, R g = 8 x 10 6 . 

FIG. 21. L2 norm convergence history of residuals for all four blocks: 

M^ = 0.5, a = 10.0°, R e = 8 x 10 6 . 

FIG. 22. Lift coefficient comparison between TNS and full-potential 
code: NACA 0012, M, = 0.5, R e =8x 10 6 . 

FIG. 23. Drag coefficient comparison between TNS and full-potential 
code: NACA 0012, M, = 0.5, R g = 8 x 10 6 . 

FIG. 24. Perspective view of particle trajectories over a stalled 
wing: NACA 0012, M^ = 0.5, a = 15.0°, R e = 8 x 10 6 (view from above the wing, 

looking downstream toward the wing leading edge). 
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FIG. 25. Particle trajectories over a stalled wing (view looking inboard 
from wing tip): NACA 0012, = 0.5, a = 15.0°, R g = 8 * 10 6 . 

FIG. 26. TNS lift coefficient comparison between subcritical (M^ = 0.5) 
and transonic (M m = 0.8) cases: NACA 0012. 

FIG. 27. Particle trajectories of the wing-tip vortex at maximum lift: 
NACA 0012, = 0.8, a = 6.0°, R g = 8 x 10^; the wing planform near the wing 

tip is shaded. 

FIG. 28. Numerical oil-flow pattern on upper wing surface: NACA 0012, 

= 0.8, o = 6.0°, R g = 8 x 10^ (N-node critical point, S-saddle critical 
point) . 

FIG. 29. Body-conforming zonal grids of a fighter aircraft configuration 
in the physical and computational spaces. 

FIG. 30. Base grid generated by an elliptic method: 60 x 20 * 24 grid 

points. 

FIG. 31. View of surface grid for the modified F-16A: (a) view from 

above; (b) view from side. 

FIG. 32. Planform view of the base grid. 

FIG. 33. A schematic view of the zonal structure for the modified F-16A. 

FIG. 34. Different wing-fuselage mesh topologies: (a) polar topology; 

(b) H-topology , 

FIG. 35. Grid clustering normal to no-slip surfaces: (a) coarse grid; 

(b) refined viscous grid. 

FIG. 36. Fine-coarse zone interface: (a) two-zone grid showing overlap 

at ABCD and EFGH planes in physical space; (b) grid point detail in the over- 
lap region in transformed space. 

FIG. 37. Pressure contours: = 0.9, a = 1.69°, R e = 4.5 x 10^. 

FIG. 38. Mach contours: = 0.9, a = 1.69°, R e = 4.5 x 10^. 
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